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Preface 

This  document  introduces  basic  refraction  error  estimation  for  an  electromagnetic  wave 
propagating  at  radio  frequencies  through  the  earth’s  atmosphere.  Appendices  contain  descriptive 
material  on  the  process,  algorithms,  and  application  of  refraction  correction  methods  to  radar 
tracking  data  at  each  of  the  participating  ranges. 

This  document  is  a  standalone  document;  however,  it  is  intended  as  a  prelude  to  the 
development  of  a  standard  for  correcting  of  radio  frequency  refraction  errors  occurring  in  radar 
track  measurements. 

This  document  was  prepared  for  the  Range  Commanders  Council  through  the  Electronic 
Trajectory  Measurements  Group  (ETMG).  For  questions  regarding  this  document,  please 
contact  the  Range  Commanders  Council  Secretariat  office. 

Secretariat,  Range  Commanders  Council 
ATTN:  CSTE-WS-RCC 
1510  Headquarters  Avenue 

White  Sands  Missile  Range,  New  Mexico  88002-5110 
Telephone:  (575)  678-1107,  DSN  258-1107 
E-mail:  usarmy.wsmr.atec.list.rcc@mail.mil 
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1.  Introduction 

As  a  prelude  to  the  development  of  a  standard  for  correcting  of  radio  frequency  (RF) 
refraction  errors,  the  ETMG  is  surveying  its  member  ranges  to  collect  information  on  the 
processes  and  algorithms  currently  used  for  correcting  these  errors. 

This  short  document  introduces  basic  refraction  error  estimation.  The  appendices  to  the 
document  contain  descriptive  material  on  the  process,  algorithms,  and  application  of  refraction 
correction  methods  to  radar  tracking  data  at  each  of  the  participating  ranges. 

2.  Refraction 

The  speed  of  propagation  of  an  electromagnetic  (EM)  wave  is  constant  only  in  a  vacuum. 
The  speed  of  propagation  through  any  other  medium  is  dependent  on  its  electrical  properties. 

The  medium  of  interest  is,  of  course,  the  earth’s  atmosphere;  however,  the  earth’s 
atmosphere  is  not  a  constant  medium:  its  properties  change  with  temperature,  pressure,  and 
humidity.  Generally,  the  speed  of  propagation  will  vary  continuously  as  the  EM  wave  travels 
from  the  radar  to  the  target  and  back  again. 

Refraction  occurs  as  the  changing  properties  of  the  atmosphere  constantly  alter  the  speed 
of  propagation.  For  a  radar  located  on  or  near  the  surface  of  the  earth  the  wave  is  continuously 
bent  downward.  The  effect  is  that  the  ray  follows  a  curved  rather  than  straight  path  through  the 
atmosphere. 1  Figure  1  demonstrates  the  effect. 


1  Physical  Science  Laboratory.  “Data  Systems  Manual,  Meteorology  and  Timing.”  Prepared  for  White  Sands 
Missile  Range  under  contract  DAAD07-76-0007,  September,  1979.  Pp31-53.  Available  on  request  to  the  RCC 
Secretariat. 
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Refraction  changes  the  apparent  altitude  of  the  object  so  that  it  seems  higher  than  it 
actually  is.  The  apparent  distance  to  the  object  is  also  affected  because  the  curved  path  is 
slightly  longer  than  a  direct  path  to  the  target. 

Refraction  correction  is  an  attempt  to  model  the  systematic  behavior  of  the  EM  wave  as  it 
moves  through  the  atmosphere  and  then  correct  its  effect  so  that  a  more  accurate  measurement  is 
obtained. 


A  starting  point  is  to  estimate  the  speed  of  the  EM  wave  as  it  travels  through  a  medium 
other  than  a  vacuum.  Basic  EM  theory2 *  shows  the  speed  of  propagation  through  a  medium 
depends  on  its  electrical  permittivity  (s)  and  magnetic  permeability  (p).  The  speed  within  the 
medium  is  then  given  by  Equation  1. 


1 

VP 


Equation  1 


The  effects  of  refraction  are  more  easily  analyzed  through  the  use  of  the  index  of 
refraction  (n).  The  index  is  defined  as  the  ratio  of  the  speed  of  propagation  of  EM  radiation  in  a 
vacuum  (c)  to  the  phase  velocity  in  the  medium  (v)  as  shown  in  Equation  2. 

c 

n  =  —  Equation  2 

v 

The  index  of  refraction  will  then  always  be  a  positive  number  greater  than  one. 
Propagation  of  a  wave  traveling  with  a  phase  velocity  near  the  vacuum  speed  of  light  will  have  a 
value  closer  to  one  than  a  wave  traveling  with  a  lower  phase  velocity. 

For  EM  waves  travelling  in  the  earth’s  atmosphere  the  index  is  generally  so  close  to  one 
that  it  is  often  easier  to  refer  to  the  index  in  parts  per  million  (ppm)  rather  than  use  the  index 
directly.  A  scaled  version  of  the  index  of  refraction,  called  refractivity,  often  appears  in  error 
correction  formulas. 

Refractivity,  designated  as  N,  is  related  to  the  index  of  refraction  by  Equation  3. 

N  =  (n—  1)'106  Equation  3 

For  example,  the  world- wide  average  refractivity  at  sea  level  is  often  given  316  N-units. 
Alternately,  the  index  itself  is  expressed  as  1.000316. 

Refractivity  is  comprised  of  two  parts:  the  contribution  due  to  atmospheric  pressure  and 
the  contribution  due  to  the  partial  pressure  of  water  vapor.  Thus,  refractivity  will  generally  be 
lower  in  dry  climates  than  in  humid  climates. 

Refractivity  decreases  as  the  height  above  sea  level  increases,  principally  because  the 
atmosphere  thins  as  altitude  increases.  For  example,  refractivity  at  4,000  feet  above  sea  level  is 
about  250-275  ppm. 

Often  the  effects  of  refraction  can  be  ignored  at  heights  above  100,000  feet  (30,480 
meters)  because  the  atmosphere  is  vanishingly  thin  at  that  height. 


2  Kraus,  John,  D.  and  Keith  R.  Carver.  Electromagnetics,  Second  Edition.  New  York:  McGraw-Hill  Publishing, 

1973.  366-367. 
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3.  Snell’s  Law 

Snell’s  Law  can  be  used  to  determine  the  change  in  the  direction  of  an  EM  wave  as  it 
crosses  a  boundary  between  dissimilar  mediums.  This  is  shown  in  Figure  2. 


The  angular  relationship  in  the  path  of  the  EM  ray  as  it  passes  from  one  medium  to 
another  is  given  by  Equation  4. 

n-tSina-L  =  n2sinai  Equation  4 

Combining  Snell’s  Law  and  the  Law  of  Sines  results  in  a  form  that  provides  a  complete 
description  of  the  ray  path  as  it  transits  through  discrete  atmospheric  layers  from  the  radar  to  a 
distant  object,  as  shown  in  Equation  5. 


n^cosOx  =  n2r2cos02  Equation  5 

This  process  can  be  extended  to  multiple  atmospheric  layers  to  provide  better  estimates 
of  the  true  path  length  and  change  in  the  direction  of  propagation.  Figure  3  illustrates  the 
geometry  for  Snell’s  Law  extended  to  multiple  stratified  layers.  Reducing  AN  leads  to  a  better 
estimate  but  requires  an  estimate  of  the  index  of  refraction  within  each  of  the  modelled 
atmospheric  layers. 
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Figure  3.  Stratified  Atmospheric  Layers 

This  simple  Flat  Earth  model  can  be  extended  to  the  spherical  model  shown  in  Figure  4. 


Figure  4.  Spherical  Earth  Model 


Here  the  atmospheric  layers  consist  of  spherical  shells  that  surround  the  earth,  extending 
from  the  surface  to  some  sufficiently  high  altitude  so  that  objects  at  heights  above  it  will 
experience  negligible  refraction. 

The  geometry  of  the  Spherical  Earth  model  includes  the  earth’s  radius  plus  site  altitude  as 
one  leg  of  a  large  triangle,  the  earth  radius  plus  site  altitude  plus  apparent  target  altitude  as 
another  leg,  and  the  observed  range  is  the  third  leg. 
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The  ideal  atmospheric  model  is  a  set  of  infinitely  thin  shells  shaped  by  the  gravitational 
field  of  the  earth,  each  accurately  reflecting  the  index  of  refraction  at  each  point  along  the  path. 

4.  Atmospheric  Model 

All  refraction  models  require  an  estimate  for  the  index  of  refraction  for  each  layer.  The 
earth’s  atmosphere  is  composed  of  a  constantly  changing  medium  resulting  in  different  indices 
of  refraction  along  the  path  of  the  EM  wave.  Bean  and  Thayer3  developed  a  standard  reference 
atmosphere  using  the  measured  surface  refractivity  and  an  exponential  decay  model  to  determine 
the  approximate  refractivity  at  any  altitude  above  the  surface,  as  shown  in  Equation  6. 

N(h)  =  Nse~k^h~h°^  Equation  6 


Where: 

N(h)  is  the  refractivity  at  height  h; 

Ns  is  the  surface  refractivity  determined  from  surface  weather  measurements; 
hO  is  the  surface  height  above  sea  level  of  the  observing  site; 
k  is  a  decay  constant  called  the  lapse  rate. 

Such  a  model  is  simplistic:  it  ignores  lateral  components  of  refraction,  turbulence, 
temperature  variations,  and  transient  inversion  layers;  however,  these  systematic  errors  are  often 
overshadowed  by  normal  variations  in  the  atmosphere,  so  they  can  often  be  safely  ignored 
(Physical  Science  Laboratory,  1979).  Corrections  for  refraction  are  usually  limited  to  just 
changes  in  path  length  (range)  and  apparent  height  (elevation). 

Model  variations  include  differences  in  the  lapse  rate,  limits  on  the  height  of  the 
atmosphere,  and  step  size  between  boundary  layers.  Typically,  the  lapse  rate  is  determined  by 
setting  the  top  of  the  atmosphere  at  some  height  above  which  refraction  is  negligible.  The  lapse 
rate  is  then  solved  such  that  it  produces  the  measured  No  at  the  surface  height  and  an  assumed 
negligible  index  of  refraction  at  the  top  of  the  atmosphere.  More  complex  models  may  use 
multiple  lapse  rates  to  more  accurately  reflect  the  different  meteorological  layers  within  the 
troposphere. 

5.  Atmospheric  Modeling  Parameters 

5.1  Earth  Model 

Refraction  correction  models  use  either  a  flat-  or  a  curved-earth  model.  Flat-earth  models 
are  appropriate  when  the  distance  to  a  tracked  object  is  such  that  the  difference  in  the  computed 
altitude  for  that  object  between  a  flat  model  and  the  curved  model  is  negligible.  The  geometry 
for  the  Flat  Earth  model  is  simpler  so  it  is  often  preferred  if  conditions  justify  its  use. 


3  Bean,  B.  R.  and  G.  D.  Thayer.  CRPL  Exponential  Reference  Atmosphere.  Washington:  U.S.  Dept,  of  Commerce, 
National  Bureau  of  Standards,  1959. 
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5.2  Lapse  Rate 

The  density  of  the  atmosphere  and  hence  its  index  of  refraction  at  a  given  height  is 
normally  modeled  as  an  exponential  decay.  There  is  a  span  of  choices  for  the  decay  constant, 
called  the  lapse  rate.  Since  the  earth’s  atmosphere  does  not  actually  follow  a  precise  exponential 
model,  meteorologists  often  use  different  lapse  rates  for  different  altitude  bands;  however,  a 
constant  lapse  rate  generally  provides  a  good  estimate  for  the  index  when  modeling  the 
atmosphere  for  RF  refraction  errors. 

One  often-used  method  of  determining  a  lapse  rate  is  to  create  a  decay  constant  that 
results  in  the  index  passing  through  the  surface  value  and  another  point  at  an  arbitrary  selected 
“top”  of  the  atmosphere.  One  often-used  choice  for  the  top  of  the  atmosphere  is  100,000  feet 
(30,480  meters)  where  the  refractivity  at  that  height  is  taken  as  3.9  ppm.  Figure  5  shows  the 
exponential  decay  in  refractivity  from  a  starting  point  near  sea  level  to  100,000  feet. 


Figure  5.  Central  Radio  Propagation  Laboratory  Exponential  Reference 

Atmosphere 


5.3  Surface  Refractivity 

Another  choice  is  the  starting  value  for  the  index  of  refraction  (or  equivalently  the 
refractivity)  at  the  surface  of  the  earth.  Standard  models  provide  a  value  at  sea  level  but  many 
instrumentation  radars  are  located  above  sea  level.  The  estimate  can  be  improved  a  little  for 
higher  elevations  if  the  refractivity-observing  site  altitude  is  used  instead  of  a  sea-level  value. 

For  example,  the  Central  Radio  Propagation  Laboratory  (CRPL)  Exponential 
Atmospheric  Reference  Model  (Bean  and  Thayer)  gives  the  sea  level  refractivity  as  313  ppm 
with  a  lapse  rate  of  0.143859km _1;  however,  when  the  measured  surface  refractivity  at  a  higher 
surface  elevation  is  found  to  be  250  ppm,  then  the  lapse  rate  is  adjusted  to  0. 125625km"1. 

Other  approaches  model  the  temperature  of  the  atmosphere  exponentially  and  then 
compute  refractivity  through  application  of  standard  pressure,  volume,  and  temperature 
formulas. 

5.4  Stratified  Boundary  Layers 

Atmospheric  density  changes  continuously  with  height  (as  well  as  with  other  less- 
dominant  factors).  For  computational  reasons  the  atmosphere  may  be  discretized  into  separate 
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layers,  each  with  a  constant  refractivity  value.  The  step  size  of  the  layer  is  an  important  choice. 
Large  step  sizes  reduce  the  computational  burden  but  they  may  not  provide  the  correct  level  of 
modeling  fidelity.  Closely  spaced  boundaries  add  computational  complexity  and  may  not 
provide  any  additional  accuracy  over  more  widely  spaced  layers. 


5.5  Water  Vapor  Contribution 

The  presence  of  water  vapor  in  the  atmosphere  is  an  important  factor  for  the  refractive 
bending  of  an  EM  wave,  especially  if  the  frequency  of  the  wave  is  in  the  RF  band.  Water  vapor 
is  measured  in  different  ways  and  this  leads  to  different  equations.  Often,  they  require  different 
physical  conversion  factors.  This  results  in  many  different  methods  for  determining  water  vapor 
contribution  to  refractivity. 

One  standard  technique  developed  by  Smith  and  Weintraub4  and  used  by  the 
International  Telemetry  Union  and  others  utilizes  the  local  atmospheric  temperature  and  the 
partial  pressure  of  water  vapor  to  compute  the  surface  refractivity.  This  is  shown  in  Equation  7. 


Ns  —  Ndry  +  Nwet  —  — 1 P  +  4810  — | 

Where: 

P  is  atmospheric  pressure  in  millibars; 

T  is  temperature  in  Kelvin; 

e  is  the  partial  pressure  of  water  vapor  in  the  atmosphere. 


Equation  7 


Here  the  constant  values  were  determined  empirically  by  Smith  and  Weintraub  based  on 
measured  dielectric  constants  for  dry  air  and  the  thermodynamic  properties  of  water  vapor. 


The  partial  pressure  of  water  vapor  represents  a  significant  contribution  to  refractivity. 
For  example,  the  refractivity  for  dry  air  at  sea  level  is  generally  taken  to  be  between  313  and  316 
ppm;  however,  a  test  range  in  Hawaii  uses  a  standard  surface  refractivity  of  359  ppm,  reflecting 
the  higher  absolute  humidity  of  its  subtropical  climate. 


6.  Ray  Tracing 

We  now  have  all  the  elements  needed  to  fully  solve  the  refraction  problem.  The 
exponential  model  of  the  atmosphere  provides  an  estimate  of  the  index  of  refraction  at  any  height 
above  an  observer.  The  Spherical  Earth  model  accounts  for  the  curved  shape  of  the  atmosphere 
shaped  by  the  earth’s  gravity.  Snell’s  Law  and  the  Law  of  Sines  allows  the  geometry  at  each 
layer  to  be  fully  described. 

We  need  only  shrink  the  atmospheric  layers  to  infinitesimal  widths  and  then  solve  the 
resulting  differential  equations.  Such  a  formulation  is  known  as  ray  tracing  because  the  solution 
traces  the  exact  path  of  the  EM  wave  as  it  propagates  to  a  distant  object. 

Bean  and  Thayer  demonstrated  how  Snell’s  Law  can  be  adapted  for  infinitesimal  changes 
in  the  index  of  refraction.  Ray  tracing  begins  by  adding  small  delta  values  to  the  Snell’s  Law 
variables. 


4  Smith,  Ernest  K.  Jr.  and  Stanley  Weintraub.  “The  Constants  in  the  Equations  for  Atmospheric  Refractive  Index  at 
Radio  Frequencies.”  Journal  of  Research  of  the  National  Bureau  of  Standards,  Vol  50,  No.  1,  January  1953. 
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Here  Equation  5  is  modified  to  include  infinitesimal  changes,  as  shown  in  Equation  8. 

n1r1cosG1  =  (ti-l  +  An )(r1  +  Ar)cos(d1  +  Ad)  Equation  8 

Note,  for  this  algorithm  Bean  and  Thayer  solved  for  the  bending  angle  x,  which  describes 
the  downward  curve  of  the  ray  at  any  point  along  its  path.  The  bending  angle  is  formed  by  the 
crossing  of  vectors  created  from  the  observed  elevation  angle  and  the  tangent  of  the  ray  at  the 
boundary  crossing.  The  elevation  correction  is  then  found  by  applying  simple  planar  geometry 
angle  relations.  A  detailed  diagram  showing  the  relationship  among  the  parameters  used  to 
develop  the  bending  equations  appears  in  Appendix  A. 

The  limit  as  the  deltas  approach  zero  yields  a  differential  equation.  If  the  equation  is 
expanded  and  the  products  of  differentials  are  eliminated  one  obtains  the  classic  expression  for 
the  bending  of  a  radio  wave,  which  is  shown  in  Equation  9. 

dn  ^ 

dr  = - cotO  Equation  9 

n 


When  integrated  over  the  path  length,  Equation  9  becomes  Equation  10. 

'n2  dn 


f'“  dn 

t1E2]  -  -  I  ~ 

Jnl  n 


COtQ 


Equation  10 


Where: 

x  is  the  bending  angle  through  which  the  ray  has  been  bent  from  its  original  direction; 
n  is  the  index  of  refraction; 

dn  is  the  infinitesimal  change  in  the  index  along  the  ray  path; 

0  is  the  observed  elevation  angle. 


Using  a  similar  approach,  the  distance  along  the  longer  curved  ray  path  is  then  found  with 
Equation  1 1 . 


fr  ndr 

RM=jro—  Equation  11 

The  equations  then  provide  for  a  direct  solution  of  the  path  of  the  bent  ray. 

6.1  Ray  Tracing  Using  Numerical  Integration 

An  attractive  procedure  is  to  compute  ray  path  angles  by  solving  the  integral  equation 
version  of  Snell’s  Law;  however,  these  integral  equations  must  be  solved  using  computationally 
intensive  numerical  techniques. 
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The  BAE  Systems5  and  NASA  Ray  Trace6  routines  use  recursive  numerical  integration, 
iterating  a  fixed  number  of  times  or  until  the  apparent  height  of  the  object  reaches  a  stationary 
value.  The  NASA  ray  trace  adds  adaptive  range  corrections  based  on  empirical  measurements. 

A  description  of  these  two  algorithms  appears  in  Appendix  B. 

6.2  Simplifying  Computational  Procedures 

Many  of  the  refraction  correction  techniques  used  at  the  test  ranges  were  developed  at  a 
time  when  the  computational  power  to  solve  complex  integral  equations  in  a  real-time 
environment  did  not  exist.  Because  the  integrals  are  computationally  difficult  to  solve  many 
refraction  models  rely  on  the  simpler  algorithms  that  provide  approximate  but  adequate 
corrections  for  the  effects  of  refraction. 

A  crude  correction  can  be  found  using  the  4/3  Earth  Radius  model.  Here  the  refractive 
index  is  assumed  to  decrease  linearly  with  height.  The  curvature  of  the  rays  is  then  such  that,  if 
plotted  in  a  geometry  where  the  earth  has  a  radius  k  times  greater  than  its  true  radius,  they  will 
appear  as  straight  lines.7  A  standard  value  for  k  is  4/3.  Such  a  model  may  provide  reasonable 
results  if  the  altitude  of  the  object  is  relatively  low.  The  4/3  approximation  is  used  only  for 
“back  of  the  envelope”  guesswork.  None  of  the  test  ranges  use  this  method  to  correct  for 
instrumentation  radar  errors. 

Where  conditions  warrant,  a  simple  Flat  Earth  model  with  fairly  wide  atmospheric  layers 
may  suffice.  A  spherical  model  must  be  used  when  the  curvature  of  the  earth  becomes 
significant  over  the  region  of  interest.  Here  the  method  is  to  compute  exact  values  for  a  few 
specific  points  in  a  range-elevation  volume  of  interest.  These  are  then  presented  in  a  coverage 
plot.  The  correction  factors  can  be  interpolated  by  eye  or  the  points  can  be  fit  to  representative 
curves  so  that  interpolations  can  be  computed  though  mathematical  algorithms. 

Work  by  NASA  in  the  mid-1960s  for  the  Gemini  Spacecraft  program  approximated  the 
solution  to  the  bending  angle  integral  with  a  simple  linear  relationship  to  the  surface  refractivity, 
the  cotangent  of  the  observed  elevation  angle,  and  a  polynomial  with  empirical  coefficients.8 
This  is  shown  in  Equation  12. 

Equation  12 

1.072014'10-2  1.279119-10-8  1.227363-10_81 

1.03585796  - + - - - 5 - 

£  £6 

Here  s  is  the  observed  elevation  angle.  The  correction  was  found  to  be  valid  for  observed 
elevations  between  2  and  10  degrees.  Above  10  degrees  only  the  first  term  in  brackets  is  used. 

Models  often  compute  a  few  salient  points  over  the  region  and  then  fit  curves  through  the 
points  to  provide  a  simple  method  to  interpolate  between  them.  Often,  empirical  corrections  are 


r  «  [(/Vs-106)(cot(V)] 


5  BAE  System,  “RF  Refraction  Correction,”  n.d.  Available  on  request  to  the  RCC  Secretariat. 

6  Agee,  B.  and  Don  Sammon.  “NASA  Ray  Trace  Algorithm,”  adapted  from  unknown  NASA  author.  Unpublished 
C  Source  Code,  September  7,  1995.  Available  on  request  to  the  RCC  Secretariat. 

7  Merrill  I.  Skolnik.  Radar  Handbook.  New  York:  McGraw-Hill,  1970,  2-45. 

8  P.  E.  Schmid.  Atmospheric  Tracking  Errors  at  S-  and  C-Band  Frequencies.  NASA  Technical  Note  TN  D-3470. 
Washington,  D.C:  National  Aeronautics  and  Space  Administration,  August  1966,  p.  9. 
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used  to  improve  the  initial  approximation.  The  simple  Kermit-Pearson  formulation  displayed  in 
Equation  13,  appears  to  be  based  on  this  approach. 


Ae  = 


kle'r'cos(el ) 
k2e  +  r  ■  sin(el ) 


Equation  13 


One  method  by  Maron9  computes  accurate  correction  values  over  a  range-elevation  grid 
and  then  fits  spline  curves  through  the  subsets  of  the  computed  points.  The  correction  for  a 
given  range-elevation  observation  pair  is  then  found  by  interpolation  between  smooth 
mathematical  curves. 

Another  approach  by  Milnarich10  is  to  recast  the  non- integral  form  of  Snell’s  Law  so  that 
it  is  suitable  for  an  iterative  solution.  The  iteration  then  computes,  in  reverse  fashion,  an 
observation  angle  error  using  an  exponential  atmosphere  model  for  the  index  of  refraction.  It 
stops  when  the  sum  of  the  geometric  and  error  correction  matches  the  actual  observed  angle. 

Variations  of  these  simplifying  techniques  lead  to  the  plethora  of  different  refraction 
correction  algorithms  used  at  the  test  ranges. 

Preliminary  comparisons  of  the  various  approaches  suggest  that  both  ray  tracing  and 
curve  fitting  approaches  can  provide  reasonably  accurate  results.  Comparisons  at  White  Sands 
Missile  Range  (WSMR)  between  the  simple  Kermit-Pearson  Flat  Earth  model  and  the 
sophisticated  NASA  Ray  trace  algorithm  suggest  that,  at  least  at  elevations  above  five  degrees, 
the  difference  between  methods  may  be  on  the  same  order  as  the  differences  caused  by  small 
variations  in  the  atmosphere.  The  general  rule  appears  to  be  that  applying  some  correction  for 
refraction  is  more  important  than  the  choice  of  the  method  used  to  implement  the  correction. 


6.3  Other  Computational  Factors 

Even  when  refraction  corrections  are  applied  to  radar  track  data,  the  test  ranges  may  use 
different  techniques  to  implement  the  corrections.  For  example,  some  algorithms  do  not  attempt 
to  estimate  refraction  when  the  elevation  is  below  zero.  Others  mirror  the  below-zero  correction 
with  its  positive  angle  counterpart  to  avoid  any  discontinuities. 

Some  algorithms  add  empirical  corrections,  especially  for  low  elevation  angles  where  the 
models  tend  to  deviate  markedly  from  the  actual  EM  propagation. 

In  addition  to  applying  refraction  to  outgoing  data,  some  ranges  apply  a  reverse  refraction 
to  incoming  cueing  data  so  that  an  instrument  will  be  cued  to  the  observed  rather  than  actual 
location  of  an  object. 

Other  procedures  call  for  applying  the  correction  only  under  certain  conditions.  For 
example,  radar  data  may  not  be  corrected  for  refraction  unless  the  radar  is  indicating  a  close-loop 
track  of  an  object. 


9  D.  E.  Maron.  “Refraction  Correction  Algorithm,”  MOTR  Program  Memorandum  MOTRPM-2000-020.  RCA 
Corp,  June  3,  1985.  Available  on  request  to  the  RCC  Secretariat. 

10  Paul  Milnarich.  “A  Method  for  Predicting  the  Apparent  Elevation  Angle  for  a  Tropospheric  Object  From  a 
Specified  Elevation  Angle  and  Height”  (dissertation,  UTEP,  1966), 
http://digitalcommons.utep.edu/dissertations/AAIEP00334/. 
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7.  Survey  of  Range  Radar  Refraction  Correction  Methods 

Appendix  C11  of  this  document  describes  the  refraction  correction  methods,  and  to  the 
extent  that  they  are  available,  the  model  parameters,  at  each  of  the  participating  test  ranges. 

8.  Future  Work 

The  long-term  goal  is  to  select  a  standard  algorithm  that  is  suitable  for  a  wide  range  of 
conditions  and  that  can  be  used  at  most  member  ranges  to  correct  radar  refraction  errors.  Future 
work  should  first  include  a  down  selection  of  suitable  candidate  algorithms,  coding  the  candidate 
algorithms  in  a  common  software  language  such  as  MATLAB,  and  then  evaluating  their 
effectiveness  on  multiple  sets  of  data  culled  from  member  ranges. 

One  candidate  algorithm  from  BAE  Systems,  the  RIR980  refraction  model,  has  already 
been  deployed  on  several  test  ranges.  Its  methods  are  well  documented  and  code  to  implement 
the  algorithm  is  readily  available.  As  such,  it  is  a  good  candidate  for  a  refraction  correction 
standard. 

A  slightly  more  sophisticated  algorithm,  the  NASA  Ray  Trace,  uses  a  similar  but  higher- 
order  integration  technique  and  it  includes  an  iterative  loop  and  empirical  adjustments  for  very- 
low-elevation  angle  observations.  A  description  of  these  algorithms  with  code  snippets  appears 
in  Appendix  B. 


11  Range  Commanders  Council.  “Member  Range  Refraction  Corrections.”  266-16  Appendix  C.  November  2016. 
Available  to  RCC  members  with  Private  Portal  access  at 

https://wsdmext.wsmr.armv.mil/site/rccpri/Publications/266-16  Survey  Radar  Refraction  Error  Corrections/266- 

16  Appendix  C.pdf. 
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APPENDIX  A 


Bending  Angle  Geometry 


r0  Earth  radius 
hs  observation  site  elevation 
h  target  altitude 
0O  observed  elevation  angle 
s  elevation  refraction  error 
P  true  elevation  angle 
Rm  measured  range 
R  ray  path 
Ro  true  range 
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APPENDIX  B 


Ray  Tracing  Techniques 


B.l  BAE  and  NASA  Ray  Trace  Algorithms 

Ranges  having  FPS-16  upgrades  developed  by  BAE  Systems  of  Fort  Walton  Beach  FL 
use  identical  refraction  correction  algorithms  (BAE  Systems).  Currently  Yuma’s  FPS-16  class 
radars  and  Eglin’s  RIR980s  utilize  this  algorithm.  Another  range,  WSMR,  is  slated  to  begin 
using  this  algorithm  in  one  upgraded  FPS-16  radar  as  soon  as  it  becomes  operational  in  January 
2017.  Four  other  upgrades  are  planned  in  the  near  future  and  all  of  them  will  use  the  BAE 
refraction  correction  algorithm. 

Eglin  reports  this  algorithm  provides  results  that  are  consistent  with  older  algorithms 
developed  for  their  FPS-16  radars.  If  other  ranges  can  confirm  its  effectiveness  then  it  could  be 
the  basis  of  an  RCC  Refraction  Correction  standard. 

The  WSMR  radar  community  also  uses  a  NASA  Ray  trace  routine  that  is  similar  to  the 
BAE  formulation.  The  exact  origins  of  this  routine  have  been  lost  but  Bill  Agee,  a 
mathematician  at  WSMR,  modified  the  original  NASA  Fortran  code  to  include  empirically 
derived  range  bending  for  targets  at  long  ranges  and  low  elevation  angles.  Don  Sammon,  also  at 
WSMR,  later  converted  the  Fortran  code  to  C. 

Both  the  BAE  and  NASA  Ray  trace  routines  use  similar  approaches  to  recursively  solve 
the  bending  angle  integral  using  a  Gaussian  Quadrature  numerical  technique.  The  BAE 
algorithm  uses  a  fourth-order  polynomial  while  the  NASA  routine  uses  an  eighth-order 
polynomial. 

The  differential  equation  for  the  bending  angle  is  shown  in  Equation  9.  Integrating  over 
the  path  length  yields  Equation  10. 

This  integral  must  be  solved  numerically.  One  could  compute  function  values  at  equally 
spaced  points  and  then  multiply  each  result  by  a  set  of  weighting  coefficients.  The  familiar 
Trapezoid  Rule  is  one  such  technique;  however,  these  ray  trace  routines  use  a  more  sophisticated 
Gaussian  Quadrature  technique  to  solve  the  integral. 

Gaussian  Quadrature  methods  provide  the  additional  degrees  of  freedom  by  selecting  not 
only  the  weighting  coefficients,  but  also  the  location  of  the  abscissas. 12  Gaussian  Quadrature 
methods  work  only  if  the  integrand  is  very  smooth  and  free  of  discontinuities. 

One  attractive  feature  of  the  Gaussian  Quadrature  is  that,  when  the  integrand  is  arranged 
as  a  polynomial  times  some  weighting  function,  and  the  weights  and  abscissas  are  carefully 
selected,  the  method  can  provide  exact  results.  Hence  the  approximation  shown  in  Equation  B-l 
can  be  exact  if  f(x)  is  a  polynomial. 


Equation  B-l 


i=l 


12  William  H.  Press,  et.  al.  “Gaussian  Quadrature  and  Orthogonal  Polynomials”  in  Numerical  Recipes  in  C. 
Cambridge  University  Press,  New  York  1988.  131-137. 
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The  selection  of  the  abscissas  at  which  the  integral  is  to  be  evaluated  corresponds  to  the 
roots  of  the  polynomial.  The  weighting  function  is  tied  to  the  theory  of  orthogonal  polynomials. 
The  determination  of  the  orthogonal  weighting  function  exceeds  the  scope  of  this  document; 
however,  numerical  values  are  tabulated  in  many  mathematical  handbooks. 

Both  routines  solve  the  integral  using  the  Gaussian  Quadrature  procedure  and  then  iterate 
several  times  to  make  sure  that  the  result  is  stable.  The  NASA  Ray  Trace  also  is  adaptive  in  the 
sense  that  it  adds  empirical  range  bending  when  the  observed  range  is  greater  than  50,000  yards 
and  the  observed  elevation  is  less  than  3  degrees. 

B.2  NASA  Ray  Trace  Algorithm 

The  C  code  for  the  NASA  Ray  Trace  program  (Agee  and  Sammon,  1995),  with 
additional  comments,  is  included  below. 

// - 

//  NASA  Ray  Trace  algorithm 

// - 

/* 

This  routine  accepts  apparent  elevation  and  range,  and 
computes  a  negative  delta  e  and  delta  r  correction 
using  a  9  point  gaussian  quadrature  integrator. 

The  algorithm  iterates  the  integral  solution  up  to  five  times 
or  until  the  change  in  the  computed  target  altitude  is  less 
than  one  foot  between  successive  interations. 

The  algorithm  also  includes  an  emperically  derived  bending 
correction  for  range  when  measured  range  is  greater  than 
50kyards  and  the  elevation  is  <3  deg. 

I/O  Parameters: 

rad  =  earth  radius  in  feet  (at  station) 
eo  =  apparent  elevation  in  radians 
ro  =  apparent  range  in  yards 

fn  =  surface  index  of  refraction  (e.g.  1.000275) 

de  =  elevation  correction  in  radians 
dr  =  range  correction  in  feet 

Algorithm  is  based  on  Bill  Agee's  (NRO-AM)  adaptation  of  the  NASA  ray 
trace  program  with  empirically  derived  range  bending  corrections 
for  elevations  below  3  degrees. 

Don  Sammon  NRO-DR-S  (575) -678-2950  7  Sep  95 

Additional  comments  added  2  June  2016 


*/ 

#def ine  pi2  1.570796327 
raytr (rad, ro, eo, fn, de, dr ) 


double 

rad; 

double 

ro ; 

double 

eo; 

double 

fn; 

double 

*de  ; 

double 

*dr ; 
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{ 

// - 

/ / declarations 

// - 

//  constants  weighting  values  for  gaussian  integration 


double 

al [9] ={ .04063719418, .09032408035, .1303053482, 

.15617353850, .16511967750, .1561735385, 

.13030534820, .09032408035, .04063719418 
} ; 

double 

u [9] ={- . 015 91 98 8 02, -.081 98444 63,- .1933142836, 

-.3378732883, -.5000000000, -. 6621267117, 

-.8066857164, -. 9180155537,-. 9840801198 
}; 

//  empircally  derived  coefficients  for  the  range  bending  correction 
double  rpl [ 1 4 ] = {  0.5, 0.7, 1.0, 1.2, 1.4, 1.6, 1.8, 


double 

2. 0,2. 5, 3. 0,4. 0,7. 0,10. ,40.0 
} ; 

po [ 14 ] =  {-.7500, -.302 8,. 1424,. 3104,. 4246,. 5058,. 56 61, 

.6127,  .6929, .7415, .7991, .8681, .8944, .9387 

} ; 

double 

pi [14]=  {  0. ,  .4134,  .4997,  .5330,  .5531,  .5648,  .5714, 

.5753, .5811, .5865, .5929, . 6017, . 6051, . 6111 
} ; 

double 

p2 [14]=  {  0.,. 019 93,. 03836,. 05136,. 0634 6,. 073 94,. 082 57, 

.08943,  .10035,  .10591,  .11160,  .11698,  .11871,  .12128 
} ; 

double 

p3 [14]=  {  .47, 0.74, 1.15, 1.36, 1.53, 1.66, 1.76, 
1.85,2.00,2.10,2.22,2.38,2.45,2.56 
} ; 

double 

double 

double 

double 

double 

double 

double 

double 

double 

int 

cse, rat, sth; 
f nml , stk; 
stn; 

a,  dn,  y2  ,  yl ,  fnj  ,  hj  ,  cl ,  c2  ,  c3 ,  temp; 
rp, ep; 

rr,poc,plc,p2c,p3c, bend; 
ste , theta, et ; 
sthp,  div¬ 
ide,  ldr;  //local  copies  of  de  and  dr 
i,  j , ib; 

//  — - 

//  entry  point  intialization  and  parameter  boundary  checking 

//  - 

//  boundary  checks:  return  if  range  is  less  than  0,if  the  elevation  is 
//  less  than  zero  or  greater  than  90  degrees,  or  index  is  <  0. 
ldr=0 . ; 
lde=0 . ; 

*dr  =  ldr; 

*de  =  lde; 

if  (eo<  0) return (1); 
if  (eo>pi2 ) return ( 1 ) / 
if  (ro<=0  ) return (1); 
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if  ( (fn  -1 .) <0 ) return  ( 1 )  ; 

//  compute  apparent  height  of  target  above  the  surface  of  the  geoid 
//  using  pythagorus  on  the  triangle  formed  by  the  earth  radaius+site 
//  altitude  as  one  leg  and  observed  range  and  elevation  as  the  second  leg. 
//  apparent  height  is  then  found  by  subtracting  the  earth  radius  from  the 
//  computed  third  leg 
rat  =  ro/rad; 

sth  =  rad* ( sqrt ( 1 . +2 . *rat*sin (eo) +rat*rat ) -1 . ) ; 

if  ( sth<0 . 1 ) return ( 1 ) ;  //  exit  if  target  altitude  is  neglibly  small 

//  compute  the  lapse  rate (exponential  decay  constant)  from  the  surface 
//  refractivity  provided  as  an  input  parameter 

//  Surface  Refractivity  is  computed  in  a  separate  procedure  using  wet  bulb, 
//  dry  bulb,  and  pressure  measurements  taken  at  the  radar  site  prior 
//  to  the  mission. 

//  The  other  value  used  to  compute  the  laspe  rate  is  an  assumed  N=3.36 
//  at  100,000  feet  altitude. 

//  The  Refractivity  at  any  altitude  is  then  N=Ns*exp ( -kh)  where  k  is  the 
//  laspe  rate  consistent  with  the  computed  Ns  and  N ( lOOkf t ) =3 . 36 

f nml  =  f n  -  1.0; 

stk  =  log(fnml/. 00000336) /100000.  ; 

// - 

/ /  iteration  loop 

//  The  algorithm  performs  up  to  five  iterations  of  the  integral 

//  ray  trace  equations  or  until  the  apparent  heght  of  the  target  changes 

//  less  than  one  foot  between  successive  iterations 

// - 

cse  =  cos (eo) ; 
for  ( j=0; j<5; j++)  { 

//  for  apparent  target  heights  up  to  1  million  feet  use  the  laspe  rate 
//  found  above  to  estimate  the  index  of  refraction  at  that  height 
//  For  heights  over  1  million  feet  index=l 
if  (sth<le6)  stn  =  1.0  +  fnml  *  exp (-sth  *  stk); 
else  stn  =  1.0; 


// - 

//  Solve  the  Ray  Trace  integral  equation  using  a 
//  9th  order  Gaussian  quadrature  integration 
//  see  Numerical  Recipes  in  C  Ch  4.5  pp  131-137 
a  =  1.0  -  stn/fn; 
dn  =  fn  -  stn; 
y2  =  0.0; 

yi  =  y2; 


f or  ( i=0 ; i<9 ; i++)  { 


fnj 

=  dn  *  u[i]  +  fnml; 

hj 

=  log (fnml/fnj ) /stk; 

cl 

=  1.0  +  a  *  u [ i ] ; 

c2 

=  1.0  +  hj/rad; 

temp 

=  al[i]/(fn  *  sqrt (cl 

yi 

=  yl  +  cl  *  c2  *  temp 

y2 

=  y2  +  temp/ cl; 

^cl*c2*c2-cse*cse) ) 


//  compute  the  true  elevation  angle  and  delta  e 

c3  =  1.0  +  sth/rad; 

ste  =  acos (fn  *  cse/ (stn*c3) ) ; 
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theta  =  ste-eo+f n*cse*a*y2 ; 

et  =  atan ( (cos (theta)  -  1 . 0/c3 ) /sin (theta) ) ; 
if  (eo>l . 57077887 )  et  =  eo; 
lde  =  et  -  eo; 

//  compute  the  delta  range  and  alpha  angle 
ldr  =  -a  *  yl  *  fn  *  fn/stk; 


// - 

//if  el<3  degrees  and  range  >  50kyards,  setup  empirical  range  bending 
if  (  (eo<0 .  05235  98775  6) && (ro>50  00  0  0 . )  )  { 
rp  =  ro  *  1 . e-6 ; 
if  (rp  >  40.0)  rp  =  4  0.0; 
ep  =  3.  -  eo*57 .2957795131; 

//  find  the  correct  range  index  to  select 

//  the  correct  coefficnets  for  the  adapative  interpolation 
f or  ( ib=l ; ib<14 ; ib++)  { 

if  (rp  <—  rpl [ ib] ) break; 

} 


//  interpolate  the  bending  using  the  selected  co-efficients 

rr  =  (rp  -  rpl [ ib-1 ] ) / (rpl [ ib]  -  rpl [ib-1]); 

poc  =  po[ib-l]  +  rr  *  (po[ib]  -  po[ib-l]); 

pic  =  pi [ib-1]  +  rr  *  (pl[ib]  -  pi [ib-1]); 

p2c  =  p2 [ ib-1 ]  +  rr  *  (p2[ib]  -  p2[ib-l]); 

p3c  =  p3 [ ib-1 ]  +  rr  *  (p3[ib]  -  p3[ib-l]); 

//  compute  bending  effect  from  the  co-efficients 
bend  =  p3c  -  exp (poc  +  pic  *  ep  +  p2c  *  ep  *  ep) ; 
bend  =  bend  *  fnml/. 00036; 
ldr  =  ldr  +  bend; 

} 

// - 

//  compute  a  new  apparent  height  above  the  geoid 
rat  =  (ro  +  ldr) /rad; 

sthp  =  rad  *  (sqrt(1.0  +  2 . 0*rat*sin (et ) +rat*rat )  -1.0); 

// - 

//  check  for  apparent  height  convergence 
dh  =  fabs(sth  -  sthp) ; 
sth  =  sthp; 
if  (dh<l .) break; 

} 

//return  the  results 
*dr=ldr ; 

* decide ; 
return ( 0 ) ; 
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APPENDIX  C 


Member  Range  Refraction  Corrections 

Appendix  C  is  available  to  RCC  members  with  Private  Page  access  at 
https://wsdmext.wsrnr.armv.mil/site/rccpri/Publications/266- 

16  Survey  Radar  Refraction  Error  Corrections/266-16  Appendix  C.pdf.  It  is  also  available 
to  US  Government  agencies  and  their  contractors  from  the  RCC  Secretariat  at  (575)  678-1107  or 
usarmy.wsmr.atec.list.rcc@mail.mil. 
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